home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / lucas_chk.cal < prev    next >
Text File  |  1995-07-17  |  13KB  |  358 lines
  1. /*
  2.  * Copyright (c) 1993 Landon Curt Noll
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * By: Landon Curt Noll
  7.  *     chongo@toad.com  -or-  ...!{pyramid,sun,uunet}!hoptoad!chongo
  8.  *
  9.  *
  10.  * primes of the form h*2^n-1 for 1<=h<200 and 1<=n<1000
  11.  *
  12.  * For all 0 <= i < prime_cnt, h_p[i]*2^n_p[i]-1 is prime.
  13.  *
  14.  * These values were taken from:
  15.  *
  16.  *    "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
  17.  *    Birkhauser, 1985, pp 384-387.
  18.  *
  19.  * This routine assumes that the file "lucas.cal" has been loaded.
  20.  *
  21.  * NOTE: There are several errors in Riesel's table that have been corrected
  22.  *     in this file:
  23.  *
  24.  *        193*2^87-1 is prime
  25.  *        193*2^97-1 is NOT prime
  26.  *        199*2^211-1 is prime
  27.  *        199*2^221-1 is NOT prime
  28.  */
  29.  
  30. static prime_cnt = 1145;    /* number of primes in the list */
  31.  
  32. /* h = prime parameters */
  33. static mat h_p[prime_cnt] = {
  34.     1, 1, 1, 1, 1, 1, 1, 1, 1, 1,            /* element 0 */
  35.     1, 1, 1, 1, 3, 3, 3, 3, 3, 3,
  36.     3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  37.     3, 3, 3, 3, 3, 3, 3, 3, 3, 5,
  38.     5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  39.     5, 5, 5, 5, 5, 5, 7, 7, 7, 7,
  40.     7, 7, 7, 7, 9, 9, 9, 9, 9, 9,
  41.     9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
  42.     9, 9, 9, 11, 11, 11, 11, 11, 11, 11,
  43.     11, 11, 11, 13, 13, 13, 13, 13, 13, 15,
  44.     15, 15, 15, 15, 15, 15, 15, 15, 15, 15,        /* 100 */
  45.     15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
  46.     15, 15, 17, 17, 17, 17, 17, 17, 17, 17,
  47.     17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
  48.     17, 17, 19, 19, 19, 19, 19, 19, 19, 19,
  49.     19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
  50.     19, 19, 21, 21, 21, 21, 21, 21, 21, 21,
  51.     21, 21, 21, 21, 21, 21, 21, 21, 23, 23,
  52.     23, 23, 23, 23, 23, 23, 23, 25, 25, 25,
  53.     25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
  54.     25, 25, 25, 27, 27, 27, 27, 27, 27, 27,        /* 200 */
  55.     27, 27, 27, 27, 27, 27, 27, 27, 27, 27,
  56.     27, 27, 27, 27, 27, 27, 27, 29, 29, 29,
  57.     29, 29, 31, 31, 31, 31, 31, 31, 31, 31,
  58.     31, 31, 31, 31, 31, 31, 31, 31, 31, 31,
  59.     33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
  60.     33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
  61.     33, 33, 33, 33, 35, 35, 35, 35, 35, 35,
  62.     35, 35, 35, 35, 35, 35, 35, 35, 35, 35,
  63.     35, 37, 39, 39, 39, 39, 39, 39, 39, 39,
  64.     39, 41, 41, 41, 41, 41, 41, 41, 41, 41,        /* 300 */
  65.     41, 41, 41, 41, 43, 43, 43, 43, 43, 45,
  66.     45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
  67.     45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
  68.     45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
  69.     45, 45, 45, 45, 45, 47, 47, 47, 47, 49,
  70.     49, 49, 49, 49, 49, 49, 49, 49, 49, 49,
  71.     49, 49, 49, 49, 49, 49, 51, 51, 51, 51,
  72.     51, 51, 51, 51, 51, 51, 51, 51, 51, 51,
  73.     51, 53, 53, 53, 53, 53, 53, 53, 53, 53,
  74.     53, 55, 55, 55, 55, 55, 55, 55, 55, 55,        /* 400 */
  75.     55, 55, 55, 55, 55, 55, 55, 55, 55, 55,
  76.     57, 57, 57, 57, 57, 57, 57, 57, 57, 57,
  77.     57, 57, 57, 57, 57, 57, 57, 57, 59, 59,
  78.     59, 59, 59, 59, 61, 61, 61, 61, 61, 61,
  79.     61, 61, 61, 61, 61, 61, 61, 61, 61, 61,
  80.     61, 63, 63, 63, 63, 63, 63, 63, 63, 63,
  81.     63, 63, 63, 63, 63, 63, 63, 63, 63, 63,
  82.     63, 63, 63, 63, 65, 65, 65, 65, 65, 65,
  83.     65, 65, 65, 65, 65, 65, 65, 65, 65, 65,
  84.     65, 65, 67, 67, 67, 67, 67, 67, 67, 67,        /* 500 */
  85.     69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
  86.     69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
  87.     69, 69, 69, 69, 69, 69, 69, 69, 69, 69,
  88.     69, 69, 71, 71, 71, 73, 73, 73, 73, 73,
  89.     73, 75, 75, 75, 75, 75, 75, 75, 75, 75,
  90.     75, 75, 75, 75, 75, 75, 75, 75, 75, 75,
  91.     75, 75, 75, 75, 75, 75, 75, 77, 77, 77,
  92.     77, 77, 77, 77, 77, 77, 77, 77, 77, 79,
  93.     79, 79, 79, 79, 79, 79, 79, 79, 79, 79,
  94.     81, 81, 81, 81, 81, 81, 81, 81, 81, 81,        /* 600 */
  95.     81, 81, 81, 83, 83, 83, 83, 83, 83, 83,
  96.     83, 83, 83, 83, 83, 83, 83, 83, 83, 83,
  97.     83, 83, 83, 83, 83, 85, 85, 85, 85, 85,
  98.     85, 85, 85, 85, 87, 87, 87, 87, 87, 87,
  99.     87, 87, 87, 87, 87, 87, 87, 87, 87, 87,
  100.     87, 87, 87, 87, 87, 87, 89, 89, 89, 89,
  101.     89, 89, 89, 89, 89, 91, 91, 91, 91, 91,
  102.     91, 91, 91, 91, 91, 91, 91, 91, 91, 91,
  103.     91, 91, 91, 91, 91, 91, 91, 93, 93, 93,
  104.     93, 93, 93, 93, 93, 93, 93, 93, 93, 93,        /* 700 */
  105.     93, 93, 93, 93, 93, 95, 95, 95, 95, 95,
  106.     95, 95, 95, 95, 95, 97, 97, 97, 97, 97,
  107.     99, 99, 99, 99, 99, 99, 99, 99, 99, 99,
  108.     99, 99, 99, 99, 99, 99, 101, 101, 101, 101,
  109.     103, 103, 103, 103, 103, 103, 103, 103, 103, 103,
  110.     103, 103, 103, 105, 105, 105, 105, 105, 105, 105,
  111.     105, 105, 105, 105, 105, 105, 105, 105, 105, 105,
  112.     105, 105, 107, 107, 107, 107, 107, 107, 107, 107,
  113.     107, 107, 107, 107, 107, 107, 109, 109, 109, 109,
  114.     109, 113, 113, 113, 113, 113, 113, 113, 113, 113,    /* 800 */
  115.     113, 115, 115, 115, 115, 115, 115, 115, 115, 115,
  116.     115, 115, 115, 115, 115, 115, 115, 119, 119, 119,
  117.     119, 119, 119, 119, 119, 121, 121, 121, 121, 121,
  118.     121, 121, 121, 121, 121, 121, 121, 125, 125, 125,
  119.     125, 125, 125, 127, 127, 131, 131, 131, 131, 131,
  120.     131, 131, 131, 131, 131, 133, 133, 133, 133, 133,
  121.     133, 133, 133, 133, 133, 133, 133, 133, 137, 137,
  122.     137, 137, 139, 139, 139, 139, 139, 139, 139, 139,
  123.     139, 139, 139, 139, 139, 139, 139, 139, 139, 139,
  124.     139, 139, 139, 139, 139, 139, 139, 139, 139, 143,    /* 900 */
  125.     143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
  126.     143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
  127.     143, 143, 143, 145, 145, 145, 145, 145, 145, 145,
  128.     145, 145, 145, 145, 149, 149, 149, 149, 149, 149,
  129.     149, 151, 151, 151, 155, 155, 155, 155, 155, 155,
  130.     155, 155, 155, 155, 155, 155, 157, 157, 157, 157,
  131.     157, 157, 157, 157, 157, 161, 161, 161, 161, 161,
  132.     161, 161, 161, 161, 161, 161, 161, 161, 161, 161,
  133.     163, 163, 163, 163, 167, 167, 167, 167, 167, 167,
  134.     167, 167, 167, 167, 167, 167, 169, 169, 169, 169,    /* 1000 */
  135.     169, 169, 169, 169, 169, 169, 169, 169, 173, 173,
  136.     173, 173, 173, 173, 173, 173, 173, 173, 173, 173,
  137.     173, 173, 173, 173, 175, 175, 175, 175, 175, 175,
  138.     175, 175, 175, 175, 175, 175, 175, 175, 175, 175,
  139.     179, 179, 179, 181, 181, 181, 181, 181, 181, 181,
  140.     181, 181, 181, 181, 181, 181, 181, 181, 181, 181,
  141.     181, 181, 181, 181, 181, 181, 181, 181, 185, 185,
  142.     185, 185, 185, 185, 185, 185, 185, 185, 187, 187,
  143.     187, 187, 187, 191, 193, 193, 193, 193, 193, 193,
  144.     193, 197, 197, 197, 197, 197, 197, 197, 197, 197,    /* 1100 */
  145.     197, 197, 197, 197, 197, 197, 197, 197, 197, 199,
  146.     199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
  147.     199, 199, 199, 199, 199, 199, 199, 199, 199, 199,
  148.     199, 199, 199, 199, 199
  149. };
  150.  
  151.  
  152. /* n (exponent) prime parameters */
  153. static mat n_p[prime_cnt] = {
  154.     2, 3, 5, 7, 13, 17, 19, 31, 61, 89,        /* element 0 */
  155.     107, 127, 521, 607, 1, 2, 3, 4, 6, 7,
  156.     11, 18, 34, 38, 43, 55, 64, 76, 94, 103,
  157.     143, 206, 216, 306, 324, 391, 458, 470, 827, 2,
  158.     4, 8, 10, 12, 14, 18, 32, 48, 54, 72,
  159.     148, 184, 248, 270, 274, 420, 1, 5, 9, 17,
  160.     21, 29, 45, 177, 1, 3, 7, 13, 15, 21,
  161.     43, 63, 99, 109, 159, 211, 309, 343, 415, 469,
  162.     781, 871, 939, 2, 26, 50, 54, 126, 134, 246,
  163.     354, 362, 950, 3, 7, 23, 287, 291, 795, 1,
  164.     2, 4, 5, 10, 14, 17, 31, 41, 73, 80,        /* 100 */
  165.     82, 116, 125, 145, 157, 172, 202, 224, 266, 289,
  166.     293, 463, 2, 4, 6, 16, 20, 36, 54, 60,
  167.     96, 124, 150, 252, 356, 460, 612, 654, 664, 698,
  168.     702, 972, 1, 3, 5, 21, 41, 49, 89, 133,
  169.     141, 165, 189, 293, 305, 395, 651, 665, 771, 801,
  170.     923, 953, 1, 2, 3, 7, 10, 13, 18, 27,
  171.     37, 51, 74, 157, 271, 458, 530, 891, 4, 6,
  172.     12, 46, 72, 244, 264, 544, 888, 3, 9, 11,
  173.     17, 23, 35, 39, 75, 105, 107, 155, 215, 335,
  174.     635, 651, 687, 1, 2, 4, 5, 8, 10, 14,        /* 200 */
  175.     28, 37, 38, 70, 121, 122, 160, 170, 253, 329,
  176.     362, 454, 485, 500, 574, 892, 962, 4, 16, 76,
  177.     148, 184, 1, 5, 7, 11, 13, 23, 33, 35,
  178.     37, 47, 115, 205, 235, 271, 409, 739, 837, 887,
  179.     2, 3, 6, 8, 10, 22, 35, 42, 43, 46,
  180.     56, 91, 102, 106, 142, 190, 208, 266, 330, 360,
  181.     382, 462, 503, 815, 2, 6, 10, 20, 44, 114,
  182.     146, 156, 174, 260, 306, 380, 654, 686, 702, 814,
  183.     906, 1, 3, 24, 105, 153, 188, 605, 795, 813,
  184.     839, 2, 10, 14, 18, 50, 114, 122, 294, 362,    /* 300 */
  185.     554, 582, 638, 758, 7, 31, 67, 251, 767, 1,
  186.     2, 3, 4, 5, 6, 8, 9, 14, 15, 16,
  187.     22, 28, 29, 36, 37, 54, 59, 85, 93, 117,
  188.     119, 161, 189, 193, 256, 308, 322, 327, 411, 466,
  189.     577, 591, 902, 928, 946, 4, 14, 70, 78, 1,
  190.     5, 7, 9, 13, 15, 29, 33, 39, 55, 81,
  191.     95, 205, 279, 581, 807, 813, 1, 9, 10, 19,
  192.     22, 57, 69, 97, 141, 169, 171, 195, 238, 735,
  193.     885, 2, 6, 8, 42, 50, 62, 362, 488, 642,
  194.     846, 1, 3, 5, 7, 15, 33, 41, 57, 69,        /* 400 */
  195.     75, 77, 131, 133, 153, 247, 305, 351, 409, 471,
  196.     1, 2, 4, 5, 8, 10, 20, 22, 25, 26,
  197.     32, 44, 62, 77, 158, 317, 500, 713, 12, 16,
  198.     72, 160, 256, 916, 3, 5, 9, 13, 17, 19,
  199.     25, 39, 63, 67, 75, 119, 147, 225, 419, 715,
  200.     895, 2, 3, 8, 11, 14, 16, 28, 32, 39,
  201.     66, 68, 91, 98, 116, 126, 164, 191, 298, 323,
  202.     443, 714, 758, 759, 4, 6, 12, 22, 28, 52,
  203.     78, 94, 124, 162, 174, 192, 204, 304, 376, 808,
  204.     930, 972, 5, 9, 21, 45, 65, 77, 273, 677,    /* 500 */
  205.     1, 4, 5, 7, 9, 11, 13, 17, 19, 23,
  206.     29, 37, 49, 61, 79, 99, 121, 133, 141, 164,
  207.     173, 181, 185, 193, 233, 299, 313, 351, 377, 540,
  208.     569, 909, 2, 14, 410, 7, 11, 19, 71, 79,
  209.     131, 1, 3, 5, 6, 18, 19, 20, 22, 28,
  210.     29, 39, 43, 49, 75, 85, 92, 111, 126, 136,
  211.     159, 162, 237, 349, 381, 767, 969, 2, 4, 14,
  212.     26, 58, 60, 64, 100, 122, 212, 566, 638, 1,
  213.     3, 7, 15, 43, 57, 61, 75, 145, 217, 247,
  214.     3, 5, 11, 17, 21, 27, 81, 101, 107, 327,    /* 600 */
  215.     383, 387, 941, 2, 4, 8, 10, 14, 18, 22,
  216.     24, 26, 28, 36, 42, 58, 64, 78, 158, 198,
  217.     206, 424, 550, 676, 904, 5, 11, 71, 113, 115,
  218.     355, 473, 563, 883, 1, 2, 8, 9, 10, 12,
  219.     22, 29, 32, 50, 57, 69, 81, 122, 138, 200,
  220.     296, 514, 656, 682, 778, 881, 4, 8, 12, 24,
  221.     48, 52, 64, 84, 96, 1, 3, 9, 13, 15,
  222.     17, 19, 23, 47, 57, 67, 73, 77, 81, 83,
  223.     191, 301, 321, 435, 867, 869, 917, 3, 4, 7,
  224.     10, 15, 18, 19, 24, 27, 39, 60, 84, 111,    /* 700 */
  225.     171, 192, 222, 639, 954, 2, 6, 26, 32, 66,
  226.     128, 170, 288, 320, 470, 1, 9, 45, 177, 585,
  227.     1, 4, 5, 7, 8, 11, 19, 25, 28, 35,
  228.     65, 79, 212, 271, 361, 461, 10, 18, 54, 70,
  229.     3, 7, 11, 19, 63, 75, 95, 127, 155, 163,
  230.     171, 283, 563, 2, 3, 5, 6, 8, 9, 25,
  231.     32, 65, 113, 119, 155, 177, 299, 335, 426, 462,
  232.     617, 896, 10, 12, 18, 24, 28, 40, 90, 132,
  233.     214, 238, 322, 532, 858, 940, 9, 149, 177, 419,
  234.     617, 8, 14, 74, 80, 274, 334, 590, 608, 614,    /* 800 */
  235.     650, 1, 3, 11, 13, 19, 21, 31, 49, 59,
  236.     69, 73, 115, 129, 397, 623, 769, 12, 16, 52,
  237.     160, 192, 216, 376, 436, 1, 3, 21, 27, 37,
  238.     43, 91, 117, 141, 163, 373, 421, 2, 4, 44,
  239.     182, 496, 904, 25, 113, 2, 14, 34, 38, 42,
  240.     78, 90, 178, 778, 974, 3, 11, 15, 19, 31,
  241.     59, 75, 103, 163, 235, 375, 615, 767, 2, 18,
  242.     38, 62, 1, 5, 7, 9, 15, 19, 21, 35,
  243.     37, 39, 41, 49, 69, 111, 115, 141, 159, 181,
  244.     201, 217, 487, 567, 677, 765, 811, 841, 917, 2,    /* 900 */
  245.     4, 6, 8, 12, 18, 26, 32, 34, 36, 42,
  246.     60, 78, 82, 84, 88, 154, 174, 208, 256, 366,
  247.     448, 478, 746, 5, 13, 15, 31, 77, 151, 181,
  248.     245, 445, 447, 883, 4, 16, 48, 60, 240, 256,
  249.     304, 5, 221, 641, 2, 8, 14, 16, 44, 46,
  250.     82, 172, 196, 254, 556, 806, 1, 5, 33, 121,
  251.     125, 305, 445, 473, 513, 2, 6, 18, 22, 34,
  252.     54, 98, 122, 146, 222, 306, 422, 654, 682, 862,
  253.     3, 31, 63, 303, 4, 6, 8, 10, 16, 32,
  254.     38, 42, 52, 456, 576, 668, 1, 5, 11, 17,    /* 1000 */
  255.     67, 137, 157, 203, 209, 227, 263, 917, 2, 4,
  256.     6, 16, 32, 50, 76, 80, 96, 104, 162, 212,
  257.     230, 260, 480, 612, 1, 3, 9, 21, 23, 41,
  258.     47, 57, 69, 83, 193, 249, 291, 421, 433, 997,
  259.     8, 68, 108, 3, 5, 7, 9, 11, 17, 23,
  260.     31, 35, 43, 47, 83, 85, 99, 101, 195, 267,
  261.     281, 363, 391, 519, 623, 653, 673, 701, 2, 6,
  262.     10, 18, 26, 40, 46, 78, 230, 542, 1, 17,
  263.     21, 53, 253, 226, 3, 15, 27, 63, 87, 135,
  264.     543, 2, 16, 20, 22, 40, 82, 112, 178, 230,    /* 1100 */
  265.     302, 304, 328, 374, 442, 472, 500, 580, 694, 1,
  266.     5, 7, 15, 19, 23, 25, 27, 43, 65, 99,
  267.     125, 141, 165, 201, 211, 331, 369, 389, 445, 461,
  268.     463, 467, 513, 583, 835
  269. };
  270.  
  271.  
  272. global    lib_debug;    /* from lucas.cal */
  273.  
  274.  
  275. /*
  276.  * lucas_chk - check the lucas function on known primes
  277.  *
  278.  * This function tests entries in the above h_p, n_p table
  279.  * when n_p is below a given limit.
  280.  *
  281.  * input:
  282.  *    high_n    skip tests on n_p[i] > high_n
  283.  *
  284.  * returns:
  285.  *    1    all is ok
  286.  *    0    something went wrong
  287.  */
  288. define
  289. lucas_chk(high_n)
  290. {
  291.     local i;    /* index */
  292.     local result;    /* 0 => non-prime, 1 => prime, -1 => bad test */
  293.     local error;    /* number of errors and bad tests found */
  294.  
  295.     /*
  296.      * firewall
  297.      */
  298.     if (!isint(high_n)) {
  299.         ldebug("test_lucas", "high_n is non-int");
  300.         quit "FATAL: bad args: high_n must be an integer";
  301.     }
  302.  
  303.     /*
  304.      * scan thru the above prime table
  305.      */
  306.     error = 0;
  307.     for (i=0; i < prime_cnt; ++i) {
  308.  
  309.         /* skip primes where h>=2^n */
  310.         if (highbit(h_p[i]) >= n_p[i]) {
  311.             if (lib_debug > 0) {
  312.                 print "h>=2^n skip:", h_p[i]:"*2^":n_p[i]:"-1";
  313.             }
  314.             continue;
  315.         }
  316.  
  317.         /* test the prime if it is small enough */
  318.         if (n_p[i] <= high_n) {
  319.  
  320.             /* test the table value */
  321.             result = lucas(h_p[i], n_p[i]);
  322.  
  323.             /* report the test */
  324.             if (result == 0) {
  325.                 print "ERROR, bad primality test of",\
  326.                     h_p[i]:"*2^":n_p[i]:"-1";
  327.                 ++error;
  328.             } else if (result == 1) {
  329.                 print h_p[i]:"*2^":n_p[i]:"-1 is prime";
  330.             } else if (result == -1) {
  331.                 print "ERROR, failed to compute v(1) for",\
  332.                     h_p[i]:"*2^":n_p[i]:"-1";
  333.                 ++error;
  334.             } else {
  335.                 print "ERROR, bogus return value:", result;
  336.                 ++error;
  337.             }
  338.         }
  339.     }
  340.  
  341.     /* return the full status */
  342.     if (error == 0) {
  343.         print "lucas_chk(":high_n:") passed";
  344.         return 1;
  345.     } else if (error == 1) {
  346.         print "lucas_chk(":high_n:") failed", error, "test";
  347.         return 0;
  348.     } else {
  349.         print "lucas_chk(":high_n:") failed", error, "tests";
  350.         return 0;
  351.     }
  352. }
  353.  
  354. global lib_debug;
  355. if (lib_debug >= 0) {
  356.     print "lucas_chk(high_n) defined";
  357. }
  358.